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Abstract 

The state-space representation of an aerodynamic vortex lattice model is considered from 
a classical and system identification perspective. Using an aerodynamic vortex model as 
a numerical simulator of a wing tunnel experiment, both full state and limited state data 
or measurements are considered. Two possible approaches for system identification are 
presented and modal controllability and observability are also considered. The theory 
then is applied to the system identification of a flow over an aerodynamic delta wing and 
typical results are presented. 
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Nomenclature 


r vortex strength 

A state matrix 

B input influence matrix 

C output influence matrix 

D direct transmission matrix 

U input matrix 

A- discrete-time eigenvalues 

A diagonal matrix of eigenvalues 

l| /. eigenvectors corresponding to eigenvalues, \ 

*¥ matrix of eigenvectors corresponding to eigenvalue matrix, A 
Y ( Markov parameters 

H(i) Hankel matrices 

X p controllability matrix 

O p observability matrix 

w downwash velocity 

8 percent error 

b m controllability vector 
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Introduction 


In recent years, several studies have been conducted at Duke University of an 
aeroelastic delta wing system. In particular, reduced order unsteady aerodynamic models 
have been used to predict the flutter, limit cycle oscillations (LCO) and gust forced 
response of such a wing. For the case of incompressible, inviscid and irrotational flow, an 
unsteady vortex lattice method was used for time domain analyses. Moreover, once such 
a model was created, it was used to create a reduced order model (ROM) which allowed a 
reduction in the order of the vortex lattice code from a thousand degrees of freedom or 
more to the order of ten degrees of freedom while retaining essentially the same accuracy 
for the representation of fluid forces acting on a wing. 

The development of such a vortex lattice model and a reduced order aerodynamic 
model based upon aerodynamic eigenmodes can be found in Ref. [1], However, the 
representation of unsteady aerodynamic flow fields in terms of global aerodynamic 
modes can be developed in a variety of ways. In Ref. [2], the authors used Proper 
Orthogonal Decomposition modes and also a system identification model for a delta wing 
to obtain a reduced order model. The current work studies the system identification of a 
vortex lattice model in greater depth with a view to developing a methodology that can be 
used with wind tunnel experimental data. The vortex lattice model is used here as a test 
bed for a numerical experiment in a continuation of the work begun in Ref. [2] . 


Delta Wing Vortex Model 

The flow about a cantilevered half-span delta wing is assumed to be 
incompressible, inviscid and irrotational. An unsteady vortex lattice method can be used 
to model this flow. A typical planar vortex lattice mesh for the three-dimensional flow 
over a wing is shown in Figure 1 . The delta wing and the wake are divided into a number 
of vortex elements. Point vortices are placed on the wing and in the wake at the quarter 
chord of the elements. At the three-quarter chord of each plate element a collocation 
point is placed for the downwash and the velocity induced by the discrete vortices is set 
equal to the downwash (fluid vertical velocity on the wing) imposed by the prescribed 
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delta wing motion or the gust field. Unsteady vorticity is shed into the wake and 
convected with the freestream velocity. The equations expressing these relationships are 
well described in Ref [1]. 


State-Space Representation 

The relationship between a Vortex Lattice (VL) model and a state-space model 
developed by system identification methods will be established in this section. For 
simplicity, we begin with the assumption that full state information is given. This allows 
us to develop a formulation to gain basic insight before pursuing further development of 
the case where only partial state data are available. 


Full-State Data or Measurement 

The VL model can be described by 1 

AT(k + 1) + BT(k) = Du(k+l) (1) 

where TfU + l) is the strength vector of the vortex at the time step /f+1, A and B are 

aerodynamic coefficient matrices, and D is a transfer matrix for determining the 
relationship between the global vortex lattice mesh and the local vortex lattice mesh on 

the delta wing itself, and u is the downwash vector. Expressions for A, B and D are 
given in Ref. [1], Assume that there is a state-space model in the form 

x(k + 1) = Ax(k) + Bu(k) 

T(k) = Cx(k) + Du(k ) ^ 

where A is the state matrix, B is the input influence matrix, C is the output influence 
matrix, and D is the direct transmission matrix. Equation (2) gives the same map 
(relationship) as shown in Eq. (1) from u to T . Is this assumption valid? 

To answer this question, let us proceed as follows. Assuming that all states are 
measurable, the output matrix C in Eq. (2) is square and invertible. Rewrite the bottom 
equation of Eq. (2) to become 
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x(k) = C- l T(k)-C- l Du{k) 


(3) 


Substituting Eq. (3) into the top equation of Eq. (2) yields 

CT T(fc + 1) = AC~ x T(k) + [B - AC~ x D]u(k) + C~ l Du(k + 1) (4) 

Premultiplying Eq. (4) by C produces 

r(£ + l) = CAC~ x T(k) + [CB - CAC- x D]u(k)+Du(k+ 1) (5) 

Comparison of Eqs. (1) and (5) results in the following equalities 
CAC~ X = -A~ l B 

D = A~ X D (6) 

CB = CAC~ X D = -A~ X BA~ X D 

From Eq. (6), it is clear that the state-space model is not unique in the sense that we have 
the freedom to choose C. Regardless of the choice of C, however, the eigenvalues of the 

VL model, i.e., eigenvalues of -A~'B , are identical to those of A (state matrix), as long 
as C is invertible. For the case where C is an identity matrix, i.e., C = / , Eq. (6) becomes 

A = -A~ l B 

D = A~'D (7) 

B = AD = -A~ x BAT'D 


Thus, we have shown that the VL model, Eq. (1), has a state-space representation, Eq. 
(2), with system matrices. A, B, and D, uniquely determined by Eq. (7), assuming that the 
output matrix C is an identity matrix. This state-space model may be used for model 
reduction and system identification. 

Assume that we are given a (discrete time) sequence of 


r = [r(i) 

r(2) • 

•• no] 

(8) 

U = [w(l) 

w( 2) •• 

• «w] 

(9) 


This may be a sequence of experimental data obtained from measurements or, as in the 
present case, numerical data from the vortex lattice model. One may use these sequences 
T and U to identify a set of system matrices, A, B, C, and D. If the identified C is an 
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invertible square matrix, the eigenvalues of the state matrix A are, indeed, the eigenvalues 

of the VL model, i.e., the eigenvalues of -A~ l B . To make C square and invertible 
requires that the number of measurements (rows) is identical to the number of states 
(colunuis). Furthermore, all rows of U must be linearly independent for the system 
identification to be valid. If the rows of U are not linearly independent, the matrix B 
cannot be properly identified. Does this mean that the eigenvalues of A may not be 
identified properly? The answer is "no". The eigenvalues of A can be fully identified if 
some of the rows can sufficiently excite the system eigenvalues. See, for example, Ref. 
[3]. 


Limited State Data or Measurements 

For simplicity, we have first developed the basic formulation based on a set of full 
state data or measurements. In practice, the spacial or temporal measurements may be 
extremely limited probably tens of spacial locations and/or time steps. For example, what 
if one is given only the first few elements of F (k) for k = 1,2,* generated by Eq. (1)? 

Can we use this very limited number of elements to identify the eigenvalues of - A~ l B ? 
One might intuitively conclude that it is very difficult, if not impossible. But, in fact, 
some progress can be mode as follows. 

Let Eq. (1) be transformed to a new form in terms of modal coordinates, i.e., 

(£ + 1) + Ar m (k) = D m u(k + 1) (10) 

where 


A = *P“ 


A~ l B 


'P = diag[A, 1 X 2 - X n ] 


r m (k) = A>- l T(k) 
D m = ¥ ~ l A~ l D 


The quantities X x X 2 ■■■ X n are the eigenvalues of the VL model, 
and')/ \\f 2 ■■■ \|/ B are the corresponding eigenvectors that form the matrix ML 
Assume that we have only one input, i.e., u(k+l) in Eq. (10) is a scalar. All the modes 
will be excited if the colunui vector D m has no zero elements in it, and n(h+l) is 
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nonzero. If some elements of D m are zero, the corresponding eigenvectors will not be 


excited. Other inputs may be added to excite such modes, of course. The elements of D m 


with relatively larger magnitude will excite the corresponding eigenmode to have a larger 
response. Assume that we have excited all the modes by properly choosing some input 
locations and their excitation signal, i.e., none of the columns in 

r ffl = [r M (i) r.( 2 ) - r m (f)] 

is a zero vector. The output data or measurements will become 

r = [r(i) r( 2) - r(f)] = ‘P[r M (i) r.(2) ••• r M oo] 


If 'P is invertible, each row of T should include all information of the VL eigenvalues 
and eigenvectors. In theory, one should be able to extract the information of eigenvalues 
from one row of T unless there exists some repeated eigenvalues. 

Let us rewrite Eq. (1) to become 


r(Jt + 1) = AT(k) + Du(k + 1) 
y(k) = C T(A:) 

where 

A = -A~'B 
D = A~ l D 


and C is an m x n matrix with m < n . Each row of C could have a unity at one location 
and zeros at all other locations so that y(k) includes only the desired elements in T ( k ) . 
From Eq. (1 1), the following matrix equation can be easily derived, i.e.. 


y(k) 
y(k+ 1) 



c 

CA 

r(k)+ 

"0 0 0 

0 CD 0 

0 " 
0 

u(k ) 
u(k+ 1) 

_y(k + p- 1)_ 


CA p ~ l 


0 ca p ~ 2 d cad 

CD_ 

u(k+ p- 1) 


( 12 ) 


or 

y p (k) = Q p T{k) + T p u p (k ) (13) 

where 
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y(k ) u(k ) 

X^ + l) n s. u(k+ 1) 

X£ + .P _ 1)_ m(^ + /?-1) 

0 0 0 0 

0 CD 0 0 

0 C4D CD 

where p is an integer that must be chosen to make 0 p a square matrix of pmxn with 
pm=n and m is the number of outputs. There is a great chance that such an integer p does 
not exist. We will discuss this case later. The quantity Y is a pm x pr matrix with r 
being the number of inputs. 

The quantity y p (k) is a pm x 1 colunm vector stacking together the m x 1 vectors, 
y(k), y(k+l), y(k+p-l) from different time steps, where m is the number of outputs. 
Similarly, the quantity u p (k) is a pr x 1 column vector consisting of rxl vectors, u(k), 
u(k+l), u(k+p-l), where r is the number of outputs. Solving Eq. (13) for T (k) yields 
r (k)=Q- p l y p (k)-Q; l r p u p (k) (14) 

Substituting T (7c) and T(T + 1) from Eq. (14) into Eq. (11) and pre-multiplying the 
resultant equation by 0 /; produces 

T/^+i)=© p i0>/^)-©^©;V,w + V,(^ +1 ) +0 ^ +1 ) ( 15 ) 

Assrmie that the mmiber of inputs is r, i.e., u(k+l) is a rxl vector. Define Y as 

Y /; = T p with the first y zero columns replaced by 0 p D (16) 

Equation (15) becomes 

y p (k+i) = e p Ae; l y p (k)-Q p AQ; l r p u p (k)+i p u p (k+\) (17) 

For easy reference, let us call Eq. (17) the Generalized Vortex Lattice (GVL) model. This 
equation has a state-space representation similar to Eq. (2), i.e., 
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x(k + 1) = Ax(k) + B p u p ( k ) 

y p (k) = C p x(k) + D p u p (k) 

where 

a = q p aq ; 1 

B p =Q p AQ-;[T p -r p ] 
c =/ 

p pmXpm 


( 18 ) 


(19) 


The output matrix C p is chosen to be an identity matrix. Equation (18) is a state-space 
model representing the map from the input vector u(k) of pr x 1 to the output vector 
y p (k) of pm x 1 . Now what we are looking for is the map from input vector u(k) of r x 1 


to the output vector y(k) of m x 1 . Careful examination of Y as defined in Eq. (16) 
reveals that 


Y p -Y p = [© pJ D 0 pmx(p _ l)r j (20) 

where 0 x(i ,_ 1)r is an pm x (p - l)r zero matrix. Equation (18) thus reduces to 
x(k + 1) = Ax(k) + Bii(k) 

y(k) = Cx(k) + Dn(k) ( ^ 

where 

a=g p aq-; 

B = Q p D 

( 22 ) 

^ _ [^bxb ^ 

D = First m rows of Q p D 

The measurement or data output equation for y(k) is obtained by taking the first m rows 
of y p (k) shown in Eq. (18) and noting that / is an identity matrix and Y ;) is a 

special matrix, i.e., the upper right-hand entries are all zeros. The matrix A is related to A 
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by a similarity transformation such that they share the same eigenvalues. We thus 
conclude that the eigenvalues for the VL model can be obtained from the state-space 
model shown in Eq. (21). The fundamental assumption is that the matrix 0 /; must be 

invertible. The matrix©^ is commonly referred to as the observability matrix in the 
control field. If some of the eigenvalues of the VL model are not observable, then 
0 /; will not be full rank and thus not invertible. Then the identified state-space matrix A 
shown in Eq. (21) would not include the un-observable eigenvalues. The maximum rank 
of 0 /; is the number of states, n. Therefore, even though we have an oversized pmxn 

matrix 0^ with pm > n , a good system identification technique should be able to identify 

a minimum-size matrix A that would include only the observable and controllable VL 
eigenvalues and eigenvectors. 


System Identification 

Several methods may be used for system identification (see Ref. [3]). Each has its 
own disadvantages and advantages. Here two simple approaches are presented, (1) 
Generalized Vortex Lattice (GVL) model identification and (2) Eigensystem Realization 
Algorithm (ERA). 


GVL Model identification 

The GVL model identification begins with Eq. (17). Let us rewrite Eq. (17) to 
become 

y P (* + !) = A y P (k) ~ A ^p u P ( k ) + ^p u p (k + 1) (23) 

where the state matrix A is defined in Eq. (22). Let the time index k run from 1 to t . 
Equation (16) will then produce the following matrix equation. 

7/2) = AY p ( 1) - AY p U p (\) + T p U p (2) 
or 
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(24) 


"r,a)‘ 

U p ( 1) 

U p (2)_ 

where Y p ( 1), Y p ( 2), U p (\) , and U p ( 2) are defined by 

y p (^) = [3 ; p(^) ^(^+1) ••• y p (t-p+k- 1)] 

X*) X^ + l) ••• y(£-p+k- 1) 

y(£+l) X^+2) ••• y(£-p + k) 

y(k+p-l ) y(A: + /?) ••• y(^+^ - 2) 

and 

= [“„(*) «,(*+!) ••• u p (£ + p+k- 1)] 

w(£) n(A:+l) ••• u(£- p+ k- 1) 

z*(£ + l) u(k + 2) ■■■ u(£-p + k ) 

w(£ + /?-l) u(k + p) ■■■ u(£ + k- 2) 

with k = 1 and 2. Note that .4 is an n x n matrix, JY is an nx pr matrix, and Y is 
also an nx pr matrix. Correspondingly, Y (1) and Y ( 2) are both pmx(£-p ) matrices, 
and U (l) and 1/(2) are both prx(£- p) matrices. 

In Eq. (24), we have three matrices. A, AY , and Y , with a total of n x (n + 2 pr) 

unknowns to be determined. The matrix U ( 1) is of the size (n + 2pr)x(£- p) . If 

U p (2)_ 

the data length £ is given such that £ - p > n + 2pr , then the maximum rank of the 
matrix is n + 2 pr if all the rows are linearly independent. Unfortunately, the matrices 
U p (1) and U (2) have all rows in common except the first row of U p {\) and the last row 

of U p ( 2). This means that the three matrices, A, AY p , and Y p cannot be solved 
uniquely. To solve the state matrix uniquely, let us define 




Y p (2) = [a -AY p Y p 
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(25) 


Y = [o„ xr -Ar p ]+\r p o„ xr 

Using Eq. (25), Eq. (24) can be rearranged to become 

r W 1 

r ' (2 n j v \uU4 (26) 

The matrix U p+l ( 1) is the union of U p ( 1) and U (2), i.e., the U ( 1) plus the last r row of 
U ( 2). Now, assume that an input signal is generated such that all rows of U p+l ( 1) are 


linearly independent, and the resulting Y (1) from the measurement also has linearly 


independent rows. Equation (26) can then be solved to yield 


r ~ - 


[ r.(i) 1 


= Y P ( 2) 

p v 7 

y p+l (i)_ 


(27) 


where f means pseudo-inverse. 


This is a least-squares solution if 


Y p ( 1) 
u p+l d) 


has more 


columns than rows. 

Equation (27) provides the state matrix A as shown in Eq. (22) which in turns 
produces the eigenvalues of the VL model. Furthermore, Eq. (27) also produces the input 

matrix B shown in Eq. (22) that is the first r columns ofY^ [see Eqs. (16) and (22)]. 


Finally, the direct transmission term D is the first m rows of B [see Eq. (22)]. 

Equation (27) for system identification is straightforward and intuitive. 
Nevertheless, it does not work for the case where the data length is so short that 


Y p ( 1) 

u p+l a) 


has more rows than columns, i.e., more spacial measurement points than time 


data steps. On the other hand, the state matrix A produced by Eq. (23) will, in general, 
not be a minimum realization in the sense that it is oversized. In other words, the integer 
p that produces the row number of Y p { 1) and 17 +1 ( 1) may be chosen larger than the 


number of states. 
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Eigensystem Realization Algorithm 

Another approach that produces a minimum-size model is based on the minimum 
realization theory. This approach uses a sequence of responses generated by a pulse input 
to the real system and then the computation and inversion of the system transfer functions 
is obtained from the input and output data (see Chapter 5 and 6 of Ref. [3]). Let 
Uj (0) = 1 (/ = l,2,...,r) and u t ( k ) = 0 (£ = 1,2,...) be substituted into Eq. (21). When the 


substitution is performed for each input element m ; (0) = 1 (i = 1 , 2 ,..., r), the results can 
be assembled into a pulse response matrix Y with dimension m by r as follows: 

Y 0 = D, Y { = CB, Y 2 = CAB, • • • , Y k = CA k ~ l B (28) 


The constant matrices in the sequence are known as Markov parameters (See Ref. [3]). 

System identification begins by forming the generalized Hankel matrix H(0), 
composed of the Markov parameters 


H( 0 ) = 


Y, 

y 2 


y 2 

y 3 


Y. 


Y, 


q + 1 


7 Y 


p + i 


Y 


p+q-l 


(29) 


The fundamental rule is that the Hankel matrix must be formed such that its rank is larger 
than the order of the system to be identified. In theory, the Hankel matrix H(0) and the 
state-space model are related by 


Y, Y 


Y 


H( 0 ) = 


Y Y 


Y 


q + 1 


Y Y 

± p 1 p + 1 


Y. 


p+q-l 


c 

CA 


[5 AB 


CA p ~ l 



(30) 


To determine A, B , C, first decompose the matrix H(0) by using singular value 
decomposition to yield 
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mp^u-Lr’ = [v r u,] S o ' ° Ik rj 

= u, I, v' + u, I if (31) 

=[i/,X“][E'T] 

where £/ and F are orthonormal matrices such that U T U = I and V' V = I , and E, is a 
diagonal matrix containing the singular values that are negligible in comparison with 
those in the diagonal matrix E r ■ Comparison of Eqs. (30) and (31) establishes the 
following equalities 

C = First m rows of \lJ r X) /:? ] 

(32) 

B = First r columns of [X) /2 F r ] 

The equalities in Eq. (32) are not unique, but they are balanced because both share the 
matrix X equally. To determine the state matrix A, another Hankel matrix must be 
formed, i.e.. 


//(1) = 


Y Y 

1 2 1 3 

Y Y 

1 3 1 A 


Y Y 

1 p + 1 ± p+l 


This matrix is formed by deleting the first column of H(0) and adding a new column at 
the end of the matrix. As a result, H(1 ) has the following relationship with the system 
matrices A, B , and C 


H(l) = 


^2 

r • 



*4 ■ 

" Y q+2 

Y P+ 1 

Y p+ 2 ■ 

" Y p+q 

C 

~ 


CA 

. r _ 



A\_B AB ••• A q - l B ] 
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Substituting Eq. (31) into Eq. (34) thus yields 

A = [uX^H( I)[x" 2 rj 

= l ; ,/2 u T r H(\)vJ i ; l/2 


( 35 ) 


The symbol f means pseudo-inverse. The orthonormal property of U and V shown in Eq. 
(3 1) has been used to derive Eq. (35). 

Assume that the state matrix A of order n has a complete set of linearly 
independent eigenvectors with corresponding eigenvalues 


which are not necessarily distinct. Define A as the diagonal matrix of eigenvalues and ¥ 
as the matrix of eigenvectors, i.e., 


A = 


\ 




X., 


(36) 


and 


^ = [¥i ¥2 ••• ¥„] 


(37) 


The identified A, B, and C can then be transformed to become A , ¥ l B , and C¥ . The 
state-space model, Eq. (21), in modal coordinates, becomes 

(*+!) = A*, (k) + B m u(k) 
y(k) = C m x m (k) + Du(k) 

where 

A = 

xjk) = K V^x(k) 

B m = ^B 

C m =cv 

The diagonal matrix A contains the modal damping rates and the damped natural 
frequencies. The matrix B m = B defines the initial modal amplitudes and the matrix 
C m = C*E the mode shapes at the sensor points. All the modal parameters of a dynamic 


system can thus be identified by the three matrices A , T 5 , and C'F . The desired 
modal damping rates and damped natural frequencies are simply the real and imaginary 
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parts of the eigenvalues A c , after transformation from the discrete-time domain to the 
continuous-time domain using the relation 

A c =-J-Log(A) (39) 


Modal Controllability and Observability 

Equation (30) shows that the Hankel matrix H(0) is formed by two matrices 
defined by 


O = 

p 


C 

CA 

CA p ~ l 


(40) 


and 

C q =[_B AB ••• A q ~ l B~\ 


(41) 


The matrix O is commonly called the observability matrix whereas C q is referred to as 


the controllability matrix. If the system possesses n states, both O and C q may have a 
maximum rank of n assuming that the integers p and q are sufficiently large. Note that 
O p and C q may not share the same rank. If either O p or C is short of rank n, then some 
of eigenvalues may not be identifiable. Let us elaborate on this statement by computing 
O p and C q in modal coordinates. 


The matrices O p and C in modal coordinates become 


r c i 


0 ¥ 

m 

C m A 




— 


C m K p -\ 


C'P[‘P“ 1 ^ v p ] p ~ 1 


and 


(42) 
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( 43 ) 


C q =[B m A B m ••• 

= 'i'-'B [ t F‘ 1 ^ l F] l p- 1 J B ••• [ 1 P- 1 ^ 1 P]‘'“ 1 

= x r l c q 

The matrix O p is the modal observability matrix whereas C q is the modal controllability 

matrix. Both matrices are coordinate dependent. Nevertheless, the Hankel matrix H(0) is 
coordinate invariant, 

W 0) = 6 p c q = O p 'V'V l C q = O p C q (44) 

There are n columns in O p and n rows in C . 

Let o j be the fth colunni of O p and cj the fth row of C q . The Hankel matrix 
H(0) can be rewritten as 

H (0) = O p C q = ^ ofif (45) 

1=1 

Both o ; and cf are associated with the ith eigenvalue . The matrix H(0) is the sum of n 
rectangular matrices formed by the product of Of] with / = Either o ; or c] will 

have contribution to H(0) only if it is not a null vector. For example, if the input actuators 
are located at the nodes of the fth mode shapes, the row vector c] is null and thus its 
norm is zero. A similar statement is also true for the output sensors and the corresponding 
vector o t . The degree of contribution of the system eigenvalue X j to H(0) is determined 

by the product of of] which yields a pmxqr matrix. The norms lloj and ||cJ| are 
commonly used to measure the degree of observability and controllability, respectively, 
for the eigenvalue L, . On the other hand, the product |o,| ||c,|| is a measure of the 

contribution of X t to H(0). A relatively larger product ||o j || |c f || means easier 

identification of the eigenvalue X i from the input and output data. Note that the input 

signal must be rich enough to excite the controllable and observable modes to be 
identified. 
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Numerical Example 


A system identification (SID) model for an unsteady aerodynamic flow has been 
created for several wing motions or gust excitations and corresponding aerodynamic 
responses. These models were derived from numerical simulations using a vortex lattice 
(VL) model for a delta wing with 55 vortex elements on the wing and 300 vortex 
elements in the wake (i.e. in Fig. 1, km = kn = 10, hnm = 20). In each case, the flow 
about the wing is excited by a certain type of prescribed downwash at the wing points, 
w(t). 

The numerical VL model produced vortex strengths at the grid points of the wing 
and in the wake, r(t), and the corresponding pressure distribution on the wing ,p(t). These 
data were then used as input for the SID model. 

The following excitations to the flow have been considered: 

1) step angle of attack, w(/)=const for t>0, 

2) frozen (fixed with respect to the fluid fixed coordinates) sharp edge gust, w(t- 
x/y)=const for (t-x/v)>0, where v is the airfoil or flow velocity, 

3) frozen gust of changing frequency, w(t-x/ y)=const sin(o) inax (t-x/ v ) 2 / 2T) , which is 
sometimes called a swept gust, where the frequency of the sweep changes from zero 
to the maximal frequency, (O _ , within the period of the sweep, T. 

4) random downwash (w at each grid point and at every time step is a random number), 

5) frozen random gust (w at the first grid point [root leading edge element] is a random 
number at every new time step, that w is then convected with the freestream velocity 
for the following panels during next steps). 

The excitation cases 1, 2, 3, and 5 are x coordinate dependent only. The excitation 
signal at a specific location x is constant along the axis y as shown in Fig. 1 . This implies 
that the excitation data matrices used for system identification are rank deficient, i.e., the 
rank is less than the number of inputs. The system matrices A. B, C, and D thus identified 
by using these types of excitation would not be accurate, in particular the colunuis in the 
input matrix B may not be independent. The number of independent columns in B is 
equal to the rank of the excitation data matrix. 
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On the other hand, the random excitation defined in case 4 is both x and y coordinate 
dependent. The excitation data matrix used for identifying the system matrices A, B, C, 
and D is of full rank (i.e., the number of inputs) with the assumption that the excitation 
data length is equal to or longer than the number of inputs. The input matrix B in 
particular would be accurately identified in theory for a controllable and observable 
system. All columns in B are linearly independent. 


Using a system identification model 

The ability of a SID model obtained from one of the excitations to predict the 
flow for another possible excitation was studied. The afonnentioned five flow excitations 
were considered. Results were obtained for either vortex strength data or pressure data. 
Consider the vortex strength data case first. Using the time histories (1000 time steps) for 
the downwash and the vortex strengths, where the latter was generated by the VL model 
for the given downwash, a SID model was obtained. That is, the state-space matrices A, 
B, C, and D for this particular downwash were computed. Then this SID model was used 
to predict the vortex strength histories due to another downwash and the results compared 
with the corresponding vortex strength time histories obtained by the VL model itself for 
that downwash. To quantify the differences between the VL and SID outputs an error, 8, 
defined as 8=100% iQ-yf/lQ / was used, where the norm / Xf is defined to be the largest 
singular value of X. In this case, both Q (VL vortex strengths) and y (SID predicted 
vortex strengths) are 1000 x 355 rectangular matrices (1000 time steps and 355 degrees 
of freedom). The same method was also used for the pressure data. There Q (VL 
pressure) and y (SID predicted pressure) are 1000 x 55 (55 elements on the wing) 
matrices. These results are presented in Figures 2 and 3 and in Tables 1 and 2. 

Consider Fig. 2. Here the percent errors are sorted based upon using the wing 
motion or gust excitation used to obtain a SID model, while the errors correspond to how 
well that SID model predicts the vortex strength (Fig 2a) and pressure (Fig. 2b) time 
histories for other excitations. A limit of 50% is used for the vertical axis scale in order to 
show small errors better. When applied to the very same downwash, the SID models 
compute the output with the less than 1% error (see Ref. [2]). As seen in the figure for 
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both the vortex strengths and pressure, the random excitation case is the best choice to 
obtain an SID model that will predict the behavior for any of the considered cases. In Fig 
2a one can see that, if not asked to predict the very special case of a random downwash, 
three other models based on sharp edge, swept gust and frozen random excitation also 
work reasonably well. 

The results in Figure 3 are sorted based upon the flow one is trying to predict with 
a SID model. For example, from the Fig. 3a one can conclude that to be able to predict 
the vortex strengths due to a frozen swept or frozen random gust, SID models obtained 
from any excitation case (except the step angle of attack excitation) would do a very good 
job, while for the prediction of vortex strengths due to a step angle of attack wing motion 
or frozen gust excitation any of considered SID models would perform well. Not 
surprisingly (after studying Fig. 2), both Fig. 3a and 3b show that only a SID model 
obtained from random downwash data can predict the time history of aerodynamic flow 
due to the random downwash. The values of these errors are presented in the Tables 1 
(vortex strengths) and 2 (pressure), where the rows of the tables relate to the Fig. 2 and 
columns to the Fig. 3. 

It is useful to consider the number of singular values kept in the SID 
models. The number following the wing motion/gust designation in the first column of 
the tables is the number of singular values retained in the identified SID model for that 
wing motion/gust. For the vortex strength data cases, the maximum number of singular 
values varies from 50 for the sharp edge gust to 300 for the random gust. This number is 
chosen such that the ratio of the largest singular value to the smallest one kept was 10 12 . 

For the vortex strength case, the total number of singular values is 300 (which, in fact, is 
equal to the total number of panels minus the number of panels on the wing). Thus, in the 
random gust case the signal is such that all the singular values are relatively high and 
close to each other magnitude. 


Lift 

As a global measure of the aerodynamic flow, consider the lift on the wing vs 
time. See Fig. 4. In Fig. 4a, b the lift time history for the motion of the wing due to the 
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step angle of attack is presented. In Fig. 4c, d the lift time history for the excitation by the 
frozen sharp edge gust is shown. (Each case is split into two figures to show in detail the 
predicted steady state lift.) The solid line is the “original” lift (obtained from the VL 
model). The left pointed and right pointed triangles represent lift time histories from SID 
models obtained from the step angle of attack motion and frozen sharp edge gust 
respectively, circles - frozen swept gust, squares - random gust, and diamonds - frozen 
random gust. (Note that the steady state lift due to either step angle of attack or sharp 
edge gust of the same amplitude is the same. That is why Fig. 4b looks very much like 
Fig. 4d.) By comparing the graphs in Fig. 4 with either results in Fig 3b or Table2, one 
can conclude that the error norm previously discussed does a good job of reflecting the 
difference between the “original” flow and those predicted by the identified models. 
However, the additional insight gained from the plots in Fig. 4 (see Fig. 4b and d) is that 
SID models from all excitations (but the swept gust) predict the steady state flow very 
well. 

From the above numerical results, it was clearly shown that the system response 
may not be reproduced using an identified model from excitations other than the random 
gust or its own excitation. Note that only the random gust could identify an accurate 
model described by the system matrices A, B, C, and D for all other inputs. Other kinds 
of excitation including the step angle, sharp edge, swept gust, and frozen random could 
only reproduce its own response, although some excitations such as swept gust may 
identify a model that gives a reasonable prediction for the responses produced by step 
angle, sharp edge, and frozen random. Nevertheless, none of the identified models 
produced by the step angle, sharp edge, swept gust, and frozen random could predict the 
response excited by the random gust. This is because the excitations other than the 
random gust did not excite all system modes in order to identify a state-space model to 
represent the system accurately. Recall that only in the case of random gust all the 
excitation (downwash) data are independent while in the four other cases only ten of 
them are independent (the downwash along the y coordinate is the same). Some modes 
can only be excited by independent downwash along the x coordinate as well as the y 
coordinate. 
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Modal Controllability and Observability 


The modal observability and controllability matrices (42) and (43) were computed 
for small numbers of outputs, 7, and inputs, U 9 respectively. The ranks of the 
controllability and observability matrices ranged between 63 and 9 1 depending on what 
output and input elements were used when building (42) and (43). For example, the rank 
of the controllability matrix in the case when there is only one input applied at the tip 
element of the wing (only the last column of the full input influence matrix, B , remains) 
was 91. Remember the system posses n = 355 states, the maximum rank can be as much 
as 355. Therefore, because the ranks of the controllability and observability matrices are 
less than 355, the system is may not be controllable and/or observable if only one input 
located at the tip element of the wing is used. Moreover, a controllability vector, b m , (see 

section 5.1.1 of Ref. [3]) was computed. This vector is fomied by the product of the 
inverse of the modal matrix with the eigenvectors of A as its column vectors and a 
column of the input influence matrix B. It was previously shown (Ref. [3]) that if b m has 

a zero element, than the control force applied at that modal node cannot control this 
mode. This means that the control force is acting exactly at the node of the mode. 

It was found that the modes corresponding to the dominant branch of the 
aerodynamic flow eigenvalues (for an earlier discussion of the eigenvalues see Ref. [2]) 
are the least controllable for the input located at the tip element of the wing. The 
components of the vector b m were scaled such that the largest component is 1 . In Figure 

5, the continuous time eigenvalues of the system are shown: in Fig. 5a those eigenvalues 
that have corresponding components of the controllability vector with values larger than 
0.1 are noted by a cross, x; in Fig. 5b similar results are shown for 0.01. As one can see, 
even in the latter case the eigenvalues of the dominant branch are still not marked. 

Values of all components of the controllability vector are shown in Figure 6a. The 
mode number is plotted along the horizontal axis, while the vertical axis shows the value 
of the components of the scaled vector b m . The modes corresponding to the dominant 

branch of the eigenvalues are numbered from 1 to 25 (the first 10 of which correspond to 
the closest to zero circle), and as seen from Fig. 6a they are the least controllable. The 
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same result persists if other locations of the input force are considered or indeed if the 
force applied at more than one location including all 55 vortex lattice panels on the wing. 

An investigation of modal observability showed that the dominant modes are the 
most observable. Observability vector components are presented in Figure 6b and the 
continuous time eigenvalues that correspond to modes with scaled observability vector 
larger than certain value are presented in Figure 7. In Fig. 7a, the components larger than 
0.4 were crossed (note the two most right circles corresponding to first 1 1 modes are 
crossed) and in Fig. 7b, all of the components that are larger than 0.2 were marked. 
Unlike the controllablity case, the modes that correspond to the dominant branch are the 
most observable. 

Qualitatively the same controllability and observability results were found when 
the identified system was obtained using considered excitations. Such results for the case 
of the random gust excitation are presented in Figure 8: the modes corresponding to the 
dominant branch are the least controllable (Fig 8a) and the most observable (Fig 8b). 

When pressure data was used to obtain A, B, and C matrices, again the system 
was found to be neither controllable nor observable with only one input located at the tip 
of the delta wing, even though the scaled controllability and observability vectors look 
different in this case: these vectors for the random gust excitation are presented in Figure 
9. Note, however, that the 55 modes that appear here do not correspond to those of the 
VL model or at least such a relation was not established. 

Modal controllability provides a measure of how difficult it may be to excite the 
respective mode by the input(s). For simplicity without losing generality, let us consider 
only one input located at the tip element of the wing. A weak modal controllability 
means that the considered mode is difficult to excite by the chosen input. One may 
consider selecting an input location or several locations that would produce a relatively 
stronger modal controllability for the modes of interest. In our example shown above, it is 
seen that the dominant branch of the eigenvalues is the least controllable in comparison 
with other branches of eigenvalues regardless of the locations of the excitation inputs. 
Nevertheless, eigenvalues in the dominant branch are the most observable. In system 
identification, both controllability and observability are equally important. A 
combination of weak controllability and strong observability may be sufficient to identify 
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the eigenvalues (modes) of interest. A conventional measure of the modal identifiability 
is the product of the modal controllability and observability for the mode of interest. An 
experimentalist before designing and performing an experiment, may need to use an 
analytical state-space model to evaluate the modal identifiability for the modes of interest 
in order to choose proper locations for excitation inputs and measurement outputs. In 
addition, the excitation signal must be rich in frequency for the modes of interest to excite 
these modes properly. 


Conclusions 

A state-space representation of a theoretical vortex lattice model has been 
developed using a system identification approach. The case of limited measurement data 
has been considered. This is done in anticipation of using the proposed system 
identification method with experimental data, e.g. a system identification model might be 
obtained from pressure (sensor) measurements on a wing in a wind tunnel. Two possible 
system identification approaches have been presented. Numerical results showed the 
importance of the choice of the wing motion or gust under consideration on building 
system identification models and predicting the vortex strengths or pressure on a wing. 
Modal controllability and observability have been discussed. Numerical results for the 
chosen delta wing vortex lattice model have shown that the modes corresponding to the 
dominant branch of the eigenvalues of the flow are the least controllable and most 
observable. By least controllable one means that these modes are difficult to excite by the 
type of downwash discussed in this paper such as the step angle of attack, frozen sharp 
edge gust, frozen gust of changing frequency, random downwash, or frozen random gust. 
On the other hand, by most observable one means that the dominant modes are easy to 
observe in terms of vortex strengths or the pressure distribution on the wing. 
Identifiability of modes depends clearly on both modal controllability and observability. 
The product of modal controllability and observability is commonly used as a measure of 
modal identifiability. 
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Tab 


e 1 : Percent errors using the vortex strength data 



Step 

Angle 

Sharp 

Edge 

Swept 

Gust 

Random 

Gust 

Frozen 

Random 

Step Angle, 67 

0.41 

3.40 

246.07 

84.11 

79.09 

Sharp Edge, 50 

11.64 

0.00 

0.21 

83.90 

0.00 

Swept Gust, 233 

0.12 

0.11 

0.05 

78.04 

0.11 

Random Gust, 300 

0.00 

0.00 

0.00 

0.00 

0.00 

Frozen Random, 1 04 

13.49 

0.07 

0.25 

89.57 

0.07 


able 2: Percent errors using the pressure data 



Step 

Angle 

Sharp 

Edge 

Swept 

Gust 

Random 

Gust 

Frozen 

Random 

Step Angle, 17 

0.05 

4.25 

168 

83.8 

90.7 

Sharp Edge, 35 

20.9 

0.66 

0.37 

320 

0.65 

Swept Gust, 52 

15.3 

15.4 

2.69 

72.9 

15.3 

Random Gust, 110 

0.42 

0.42 

0.07 

0.39 

0.42 

Frozen Random,45 

25.9 

1.48 

0.74 

290 

1.46 
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Figure 1: Vortex lattice model of unsteady flow about a delta wing. 
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Figure 2: Percent error for (a) vortex strength and (b) pressure time histories for the five 
described excitations (marked by the numbers). The SID models used to simulate these 
histories are indicated along the horizontal axis. 
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Figure 3: Percent error for (a) vortex strength and (b) pressure time histories that were 
simulated by the SID models (marked by the numbers). The flows under consideration 
are indicated along the horizontal axis. 
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Figure 4: Lift time history for the excitations of (a and b) step angle of attack and (c and 
d) frozen sharp edge gust. 
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Imaginary Part of Eigenvalues Imaginary Part of Eigenvalues 



Real Part of Eigenvalues 


Figure 5: System eigenvalues. The crossed eigenvalues correspond to modes for which 
components in the scaled controllability vector are larger then (a) 0.1, (b) 0.01. The 
largest controllability vector component is 1 (cf. Fig. 6a). 
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Figure 6: Contollability (a) and observability (b) vector components. In this case, the 
state, A, and input influence, B, matrices were derived in the close form from the VL 
model with the use of Eqs. (7), where the output influence matrix, C, is an identity 
matrix. The only input is located at the tip of the delta wing. 
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Imaginary Part of Eigenvalues Imaginary Part of Eigenvalues 



Real Part of Eigenvalues 


Figure 7: System eigenvalues. The crossed eigenvalues correspond to modes for which 
components in the scaled observability vector are larger then (a) 0.4, (b) 0.2. The largest 
observability vector component is 1 (cf. Fig. 6b). 
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Figure 8: Contollability (a) and observability (b) vector components. In this case, matri- 
ces A, B, and C were identified from the vortex strength data when the wing was excited 
by the random gust. The only input is located at the tip of the delta wing. 
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Figure 9: Contollability (a) and observability (b) vector components. In this case, ma- 
trices A, B, and C were identified from the pressure strength data when the wing was 
excited by the random gust. The only input is located at the tip of the delta wing. 
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